% 石墨烯能带
clc;clear;
Gamma=[0,0];
M=[pi,pi/sqrt(3)];
K=[2*pi/3,2*sqrt(3)*pi/3];

kn=100; %每段路径k点个数
n=3; % 有几段路径

% Gamma-M段
kx(1:kn)=linspace(0,M(1),100);
ky(1:kn)=linspace(0,M(2),100);
% M-K段
kx(kn+1:2*kn)=linspace(M(1),K(1),100);
ky(kn+1:2*kn)=linspace(M(2),K(2),100);
% K-Gamma段
kx(2*kn+1:3*kn)=linspace(K(1),0,100);
ky(2*kn+1:3*kn)=linspace(K(2),0,100);

for i=1:kn
    vals=eig(hamiltonian(kx(i),ky(i)));
    E1(i)=vals(1);
    E2(i)=vals(2);
    kv=[kx(i),ky(i)];
    k(i)=norm(kv); % 第一段k路径 
    for i=kn+1:2*kn
        vals=eig(hamiltonian(kx(i),ky(i)));
        E1(i)=vals(1);
        E2(i)=vals(2);
        kv=[kx(i)-M(1),ky(i)-M(2)];
        k(i)=norm(kv)+norm(M-Gamma); % 第二段k路径
        for i=2*kn+1:3*kn
            vals=eig(hamiltonian(kx(i),ky(i)));
            E1(i)=vals(1);
            E2(i)=vals(2);
            kv=[kx(i)-K(1),ky(i)-K(2)];
            k(i)=norm(kv)+norm(M-Gamma)+norm(K-M); % 第三段k路径
        end
    end
end

%% 画图
figure();
grid on
hold on
box on % 加边框
plot(k,E1,'r')
plot(k,E2,'b')
xlim([0 k(3*kn)]) % x轴范围
xticks([0 k(kn) k(2*kn+1) k(3*kn)])
xticklabels({'\Gamma','M','K','\Gamma'})
ylabel("Energy")
title("石墨烯能带")
print("石墨烯",'-dpng','-r600') % 保存为png格式，分辨率600